Profiling of differentially expressed circRNAs and functional prediction in peripheral blood mononuclear cells from patients with rheumatoid arthritis

Abstract Background Rheumatoid arthritis (RA) is a chronic autoimmune disease associated with an increased risk of death, but its underlying mechanisms are not fully understood. Circular RNAs (circRNAs) have recently been implicated in various biological processes. The aim of this study was to investigate the key circRNAs related to RA. Methods A microarray assay was used to identify the differentially expressed circRNAs (DEcircRNAs) in peripheral blood mononuclear cells (PBMCs) from patients with RA compared to patients with osteoarthritis (OA) and healthy controls. Then, quantitative real-time PCR was applied to verify the DEcircRNAs, and correlations between the levels of DEcircRNAs and laboratory indices were analysed. We also performed extensive bioinformatic analyses including Gene Ontology (GO), Kyoto Encyclopedia of Genes and Genome (KEGG) pathway and potential circRNA–miRNA–mRNA network analyses to predict the function of these DEcircRNAs. Results A total of 35,342 and 6146 DEcircRNAs were detected in RA patients compared to controls and OA patients, respectively. Nine out of the DEcircRNAs in RA were validated by real-time PCR. There were correlations between the levels of DEcircRNAs and some of the laboratory indices. GO analyses revealed that these DEcircRNAs in RA were closely related to cellular protein metabolic processes, gene expression, the immune system, cell cycle, posttranslational protein modification and collagen formation. Functional annotation of host genes of these DEcircRNAs was implicated in several significantly enriched pathways, including metabolic pathways, ECM–receptor interaction, the PI3K–Akt signalling pathway, the AMPK signalling pathway, leukocyte transendothelial migration, platelet activation and the cAMP signalling pathway, which might be responsible for the pathophysiology of RA. Conclusions The findings of this study may help to elucidate the role of circRNAs in the specific mechanism underlying RA. Key messages Microarray assays showed that a total of 35,342 and 6146 DEcircRNAs were detected in RA patients compared to controls and OA patients, respectively. Nine out of the DEcircRNAs in RA were validated by real-time PCR, and the levels of the DEcircRNAs were related to some of the laboratory indices. GO analyses revealed that the DEcircRNAs in RA were closely related to cellular protein metabolic processes, gene expression, the immune system, etc. Functional annotation of host genes of the DEcircRNAs in RA was implicated in several significantly enriched pathways, including metabolic pathways, ECM–receptor interaction, the PI3K–Akt signalling pathway, etc.


Introduction
Rheumatoid arthritis (RA) is a chronic autoimmune disease characterized by synovial hyperplasia, pannus formation and joint damage that affects approximately 0.5-1% of adults worldwide [1]. A recent study reported that despite decreasing mortality rates due to advancements in RA management in recent years, RA continues to be linked to an increased risk of death [2]. Moreover, although major progress in the diagnosis and treatment of RA has been made in recent years, damage to articular cartilage and bone and long-term disability in RA patients are still common [3]. The pathogenesis of RA is complex and involves environmental factors that trigger disease in genetically susceptible individuals [4]. The most effective therapeutic approach of RA requires early diagnosis and adoption of an optimal nonpharmacological and pharmacological treatment, associated with periodic evaluation of therapeutic efficacy and safety [5]. Therefore, it is necessary to comprehensively elucidate the inflammatory and immunoregulatory mechanisms that underlie the etiopathogenesis of RA [6].
As a new class of endogenous noncoding RNAs, circular RNAs (circRNAs) have recently attracted more attention due to their pivotal role in the modulation of gene expression at the transcriptional and posttranscriptional levels [7]. Accumulating studies have implicated the biological functions of circRNAs such as microRNA (miRNA) sponging, regulating target gene transcription and forming RNA-protein complexes [8,9]. Importantly, several studies have also demonstrated that circRNAs may be involved in the process of mediating inflammation and immune regulation in autoimmune diseases [10][11][12]. A number of circRNAs have been shown to be aberrantly expressed in RA patients and may be involved in the pathogenesis of RA [13][14][15]. Although these lines of evidence indicate that circRNAs have the potential to serve as novel diagnostic biomarkers and therapeutic targets for RA, the expression profiles and biological functions of circRNAs in RA remain elusive.
We aimed to investigate the profiles of circRNA expression in RA in the Northwest Chinese Han population. We used a circRNA microarray assay to identify the differentially expressed circRNAs (DEcircRNAs) in peripheral blood mononuclear cells (PBMCs) from patients with RA compared to healthy controls and patients with osteoarthritis (OA). Then, quantitative real-time PCR (qRT-PCR) was applied to verify the DEcircRNAs. We next performed extensive bioinformatic analyses including Gene Ontology (GO), Kyoto Encyclopedia of Genes and Genome (KEGG) pathway and potential circRNA-miRNA-mRNA network analyses, to predict the function of DEcircRNAs, and to reveal the possible pathogenesis of RA. Our results provide a novel perspective for the pathogenesis and diagnostic markers of RA.

Study population
Patients with RA were selected from the Department of Rheumatology and Immunology of the Second Affiliated Hospital of Xi'an Jiaotong University during the period from January 2019 to December 2020. Ageand gender-matched healthy subjects and patients with OA were from the Health Examination Center and the Department of Orthopaedics of our hospital, respectively. Both healthy subjects and OA patients were recruited as the controls of RA patients. RA patients fulfilled the criteria of the 2010 American College of Rheumatology (ACR). Patients with RA and OA were excluded if they had other autoimmune diseases, haematologic diseases, malignancies or infections; had any history of other chronic diseases, such as diabetes mellitus, dyslipidaemia, thyroid dysfunction, or severe liver or kidney impairment; or received corticosteroid treatment within the last 3 months. Fresh circulating blood samples (4 mL) were collected from all of the participants, and serum was collected and stored at -80 C. This study was performed in accordance with the Declaration of Helsinki and was approved by the Research Committee of Human Investigation of Xi'an Jiaotong University Health Science Center (2018-1602), and written informed consent of the study was obtained from all participants.

Peripheral blood mononuclear cell isolation and RNA extraction
Peripheral venous blood samples were collected from all participants. Peripheral blood mononuclear cells were isolated by density gradient centrifugation using Lymphocyte Separation Medium (TBD, Tianjin, China), and were stored at -80 C in TRIzol (Invitrogen, Carlsbad, CA). Total RNA was extracted from PBMCs using TRIzol reagent and purified with a Qiagen RNA easy Mini Kit (Qiagen, Hilden, Germany) according to the manufacturer's instructions. The concentration and purity of RNA were determined using a NanoDrop ND-1000 spectrophotometer (Thermo Fisher Scientific, Waltham, MA). The integrity of RNA was assessed by electrophoresis on a denaturing agarose gel.

Microarray profiling and data analysis
The expression profiles of circRNA were detected by Agilent human circRNA Array V2.0 according to the recommended procedures provided by CapitalBio Technology Corporation (Beijing, China). The assay included the process of RNA labelling, hybridization, scanning and data analysis. Raw data were extracted and preprocessed by using Feature Extraction software V10.7 (Agilent Technologies, Santa Clara, CA). Then, normalization, quality control and differential analysis of the data were performed using Agilent Gene Spring V13.0 software. Fold-change ! 2.0 and a p value <0.05 were considered the criteria for identifying differentially expressed genes.

qRT-PCR assay validation
qRT-PCR assays were performed to validate the DEcircRNAs that were detected in the microarray analysis. Total RNA from PBMCs was extracted using TRIzol reagent (Ambion, Austin, TX) according to the manufacturer's protocol. The quantity of RNA was determined using a NanoDrop2000 (Nano Drop Products, Wilmington, DE). cDNA was synthesized using Evo M-MLV RT Premix (AG11706, ACCURATE BIO-TECHNOLOGY, Hunan, China). The qRT-PCR assay was conducted using a SYBR V R Green Premix Pro Taq HSqPCR Kit (AG11701, ACCURATE BIOTECHNOLOGY, Hunan, China) in an ABI Prism 7500 Real-Time PCR System (Applied Biosystems, Foster City, CA). The comparative CT (2 -DDCT ) method was used to obtain the fold change in circRNA expression levels. All experiments were performed in triplicate.
The interaction between miRNAs and mRNAs was predicted based on MiRanda and TargetScan (http:// www.targetscan.org). The visual circRNAs-miRNAs or circRNAs-miRNAs-mRNAs regulatory networks were established by using Cytoscape software (http://www. cytoscape.org). The function of these target genes, including biological process (BP), molecular function (MF) and cellular component (CC) were defined by GO (http://www.geneontology.org). Kyoto Encyclopedia of Genes and Genomes analysis was also performed to illuminate their involved biological pathways.

Statistical analysis
Continuous parameters were described as the means ± standard deviations and qualitative parameters as numbers (%). The differences in continuous variables were analysed using one-way analysis of variance or Student's t-test, while Chi-squared tests were performed to compare the differences in qualitative variables. The correlations between DEcircRNAs and Laboratory indices were analysed by Spearman's correlation test. All these statistical analyses were performed using SPSS software (version 16.0, Chicago, IL) with p < 0.05 considered statistically significant.

Baseline characteristics of the study subjects
This study included a total of 55 patients with RA, 55 patients with OA and 55 controls. The demographic characteristics of the subjects are given in Table 1. There were no differences with regard to age, sex or BMI among the three groups for either the sequencing cohort or the validation cohort.

Overview of circRNAs expression profiles
A total of 170,340 circRNAs were detected in this microarray analysis. Among these circRNAs, 35,342 were differentially expressed in RA patients compared with healthy controls (fold change ! 2 and p < 0.05), including 20,357 upregulated circRNAs and 14,985 downregulated circRNAs (Figure 1(A)). Additionally, 6146 circRNAs (3109 circRNAs upregulated and 3037 circRNAs downregulated) were identified as being significant in RA compared with OA using the cut-off criteria of !2.0-fold changes and p < 0.05 (Figure 1(C)). The DEcircRNAs were widely distributed across almost all human chromosomes, including the sex chromosomes ( Figure 1(A,C)). The related length distribution of these circRNAs is shown in Figure 1(B,D). Hierarchical clustering showed that circRNA expression profiles were distinctly different between RA and healthy controls (Figure 2(A)) and between RA and OA ( Figure 2(C)). Volcano plot analysis was used to visualize the DEcircRNAs (Figure 2(B,D)). The top 10 upregulated and downregulated circRNAs between RA and healthy controls and between RA and OA are displayed in Tables 2 and 3, respectively.

Functional analysis of differentially expressed circRNAs
We first performed Venn diagram analysis (Figure 3(A)) to find 3868 unique DEcircRNAs that were relevant to the pathological process of RA. Then, GO analysis and KEGG pathway enrichment were used to investigate the putative functions and cellular signalling pathways of these DEcircRNAs in RA. The GO terms are divided into three domains: BP, CC and MF. The top 30 GO processes of each domain and KEGG signalling pathways of unique and upregulated circRNAs in RA are listed in Figure 4, while those of unique and downregulated circRNAs in RA are listed in Figure 5. The results showed that the most significantly enriched KEGG pathways of these upregulated circRNAs in RA were associated with gene expression, the immune system, signalling by Rho GTPases, metabolism of proteins, signal transduction, haemostasis and adaptive Table 1. Demographic characteristics of the study subjects.

Sequencing cohorts
Validation cohorts Age (years ± SD) 53. 16    immune system. The most significantly enriched KEGG pathways of these downregulated circRNAs in RA were related to metabolism, gene expression, cell cycle, posttranslational protein modification and collagen formation.
Identification of the circRNA-miRNA-mRNA network CircRNAs have been demonstrated to play a vital role in the transcriptional regulation of genes by functioning as miRNA sponges. Therefore, we identified the miRNAs targeted by the top 10 upregulated (red nodes) and downregulated (green nodes) circRNAs in RA patients compared with healthy controls and OA patients ( Figure 6(A,B)). The potential miRNA targets of these DEcircRNAs were analysed by using the TargetScan and MiRanda databases. Then, the network diagram of the association between circRNAs and miRNAs and mRNA was generated by using Cytoscape. As shown in Figure 6(A,B), two of the upregulated circRNAs, hsa_circ_0000579_CBC1 and hsa_circ_0004417_CBC1, could target multiple miRNAs. Furthermore, we performed mutually targeted MRE enrichment (MuTaME) analysis to explore the regulatory relationships among several validated DEcircRNAs, their targeted miRNAs and the relevant mRNAs in PBMCs from RA patients. The circRNA-miRNA-mRNA   Figure 7(A,B), respectively. The regulatory function of these DEcircRNAs partially overlapped in this network analysis. Notably, RICTOR was illustrated to be the common targeted mRNA of these upregulated circRNAs. Bioinformatics analysis showed that there were binding sites between the four upregulated circRNAs and the corresponding miRNAs that could also subsequently bind to RICTOR mRNA ( Figure 8).

Discussion
Accumulating data have demonstrated that circRNAs play a critical role in regulating the immune response and inflammation [16]. RA is a chronically progressive autoimmune disease [17], and the mechanisms underlying the role of circRNAs in the pathogenesis of RA remain to be illustrated. In this study, we profiled DEcircRNAs in PBMCs from patients with RA compared to healthy controls and patients with OA. We then performed bioinformatic analyses of these DEcircRNAs and found that these circRNAs were related to metabolic pathways, ECM-receptor interactions, the PI3K-Akt signalling pathway, protein processing in the endoplasmic reticulum, leukocyte transendothelial migration, ubiquitin-mediated proteolysis and so on. We conducted qRT-PCR to verify the expression of 10 selected circRNAs in PBMCs from patients with RA. Further study of these DEcircRNAs could lead to the discovery of novel biomarkers for the early diagnosis of RA or the development of novel therapeutic strategies to prevent disease progression.
Recent advances have led to better diagnostic criteria, improved serologic testing, novel new drugs and better guidelines to manage patients with RA. However, there is no standard definition of early RA [18]. Early RA known as 'preclinical RA' clearly begins years to months before it manifests as polyarthritis. A great proportion of these patients with preclinical RA could develop increasing joint pain and swelling in the ensuing months to years, particularly in those with positive anti-citrullinated protein antibody [19]. Since early diagnosis and treatment could prevent the progression of joint damage in 90% of patients with early RA [20], it is important to identify patients with RA as soon as possible [21]. More recent studies have demonstrated that noncoding RNAs primarily including miRNA, long noncoding RNA (lncRNA) and circRNA play a critical role in inflammation and autoimmune regulation and that they are identified as promising biomarkers for the diagnosis and treatment of RA [22]. In this study, we identified DEcircRNAs in PBMCs from RA patients compared with healthy controls and OA patients. Future studies investigating the role of these DEcircRNAs in distinguishing preclinical RA from confirmed RA would be of great interest to understand the potential application of these DEcircRNAs as novel biomarkers of early diagnosis and therapeutic targets in RA patients.
In this study, functional enrichment analyses were also performed to identify the potential roles of DEcircRNAs in the pathogenesis of RA. The GO analysis of the mRNAs targeted by DEcircRNAs showed that CC organization or biogenesis, cellular protein metabolic processes, extracellular matrix component, ATP binding, adenyl ribonucleotide binding, protein binding and protein transport are implicated in the biological and cellular processes related to the pathogenesis of RA. Recently, it has been well recognized that metabolic stress is closely involved in the pathophysiology of systemic autoimmune diseases [23]. For example, glutamine metabolism was demonstrated to have distinct roles in promoting Th17 but constraining Th1 and CTL effector cell differentiation [24]. Another study found that the metabolic stress sensing protein kinase GCN2 played a critical role in regulating the tolerogenic response to apoptotic cells and limiting autoimmunity [25]. Collectively, these aberrantly expressed circRNAs in PBMCs are implicated in many pathophysiological processes of RA, especially those  regarding inflammation and immunity, and may be promising biomarkers and therapeutic targets of RA.
The KEGG pathway analysis indicated that these DEcircRNAs might be responsible for the pathophysiology of RA by regulating metabolic pathways, ECM-receptor interactions, the PI3K-Akt signalling pathway, the AMPK signalling pathway, regulation of the actin cytoskeleton, protein processing in the endoplasmic reticulum, leukocyte transendothelial migration, platelet activation, the cAMP signalling pathway and the Ras signalling pathway. All of these pathways are associated with the regulation of immune responses, and previous studies have revealed that both innate and adaptive immune dysregulation are involved in the pathogenesis of RA [26][27][28]. The AMPK cascade signalling pathway plays a central role in regulating many important pathophysiological processes, including glucose homeostasis, lipid metabolism, mitochondrial biosynthesis and protein synthesis [29], so it has become an attractive therapeutic target for many diseases. We previously found that the dopamine D3 receptor on mast cells could alleviate inflammation in mouse RA through the mTOR/AKT/AMPK signalling axis [30]. Interleukin (IL)-17 is a key regulator of immune and inflammatory responses [31], and the PI3K/Akt signalling pathway was implicated in the overproduction of IL-17 in RA patients [32]. The activation of the PI3K/Akt signalling pathway was also responsible for the neoangiogenesis of synovial tissues in patients with RA [33]. However, further study is needed to confirm this gene pathway in the pathologic process of RA.
As a new type of noncoding RNA, circRNAs function as competing endogenous RNAs of their targeted miRNAs to upregulate the expression of target genes [34,35]. In this study, we also constructed a ceRNA network to further explore the miRNA sponge role of circRNA in regulating gene expression. From the visible circRNA-miRNA-mRNA network of our results, it can be concluded that circRNAs may regulate mRNA expression by sponging one or several targeted miRNAs. For example, Rictor has been demonstrated to have an important effect on immune responses by regulating the differentiation and function of immune cells such as T cells, B cells and macrophages [11,36,37]. Therefore, Rictor plays a key role in the pathogenesis of autoimmune diseases. In this study, it was shown that hsa_circ_0054223_CBC1, hsa_circ_0053881_CBC1, hsa-circRNA13773-35_CBC1 and hsa_circ_0024203_CBC could all bind to their corresponding miRNAs to affect the expression of RICTOR. The interactive regulation among circRNAs, miRNAs and mRNAs indicated in our study could provide a novel perspective for the pathogenesis of RA.
There were some limitations in this study. First, we had relatively few RA patients for the real-time quantitative PCR (RT-qPCR) verification of the DEcircRNAs, and therefore, it is necessary to further validate our results in a larger RA cohort. We need to corroborate the association of these circRNAs with the disease activity of RA in our future study. We will also investigate the translational potential of these DEcircRNAs in the diagnosis, prognosis and therapeutic evaluation of RA patients in our future study, on account of these findings. Second, the function of the DEcircRNAs was only based on bioinformatics analysis in this study. We should perform more experiments by using methods of controlling the expression of these circRNAs to verify the precise functions of the DEcircRNAs in RA. Finally, PBMCs contain multiple cell types, including lymphocytes, neutrophils, monocytes and NK cells. Consequently, it may be worthwhile to clarify the expression profiles of circRNAs in different types of PBMCs, which deserves to be done in our future study. Given all these points, it might be necessary to investigate novel strategies for the early diagnosis, prognosis and evaluation of therapies by targeting these DEcircRNAs found in this study as laboratory biomarkers of RA and their relevant key signalling pathways in our future study.